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Intricate interplay between the periodicity of the lattice structure and that of the cyclotron motion 
gives rise to a well-known self-similar fractal structure of the energy eigenvalue, known as the 
Hofstadter butterfly, for an electron moving in lattice under magnetic field. Evolving from the 
n — Landau level, the central band of the Hofstadter butterfly is especially interesting since it 
may hold a key to the mysteries of the fractional quantum Hall effect observed in graphene. While the 
entire Hofstadter butterfly can be in principle obtained by solving Harper's equations numerically, 
the weak-field limit, most relevant for experiment, is intractable due to the fact that the size of the 
Hamiltonian matrix, that needs to be diagonalized, diverges. In this paper, we develop an effective 
Hamiltonian method that can be used to provide an accurate analytic description of the central 
Hofstadter band in the weak-field regime. One of the most important discoveries obtained in this 
work is that massless Dirac particles always exist inside the central Hofstadter band no matter how 
small the magnetic flux may become. In other words, with its bandwidth broadened by the lattice 
effect, the n = Landau level contains massless Dirac particles within itself. In fact, by carefully 
analyzing the self-similar recursive pattern of the central Hofstadter band, we conclude that massless 
Dirac particles should occur under arbitrary magnetic field. As a corollary, the central Hofstadter 
band also contains a self-similar structure of recursive Landau levels associated with such massless 
Dirac particles. To assess the experimental feasibility of observing massless Dirac particles inside 
the central Hofstadter band, we compute the width of the central Hofstadter band as a function of 
magnetic field in the weak-field regime. 



I. INTRODUCTION 

Observing the behavior of electrons in graphene un- 
der high magnetic field has played an important role not 
only for uncovering new quantum Hall states, but also 
for proving the very existence of massless Dirac parti- 
cles W, ^ . Affected by the linear dispersion near Dirac 
points. Landau levels are formed in graphene such that 
their energy is scaled as sgii{n) ^/\n\ in units of y/2hvF/lB 
with n, the Landau level index, allowed for all integers in- 
cluding positive, zero, and negative [3]. In the above, vp 
is the Fermi velocity at the Dirac point and / b = \/hc/eB 
is the magnetic length. The n = Landau level offers a 
particularly intriguing departure from the usual quantum 
Hall effect (QHE) in that its Hall coefficient is shifted by 
half an integer. With both spin and valley degeneracy 
taken into account, the consequent Hall conductance is 
predicted to be quantized in the form of 4(n -I- 1/2) in 
units of e^/ft., which exhibits beautiful agreement with 
experiment [U [2] ■ 

There is, however, a glaring omission in the discussion 
so far. In the above, the effect of lattice is completely ig- 
nored except that the electron dispersion becomes linear 
near Dirac points. The question is how valid this assump- 
tion can be. More specifically, will there be any changes 
in the Landau-level structure once the effect of lattice is 
better incorporated? Naively speaking, since the devia- 
tion from the linear dispersion occurs in relatively high 
energy, one may expect that the Landau levels should 
be more or less the same as before so that they remain 
as flat bands. In particular, the n = Landau level is 
then expected to remain as a flat band pinned exactly 
at zero energy due to the particle-hole symmetry. Seem- 



ingly innocuous, if true, this expectation gives rise to 
a very puzzling question: what determines which states 
within the n — Landau level evolves into the particle 
(or the positive energy) branch and which into the hole 
(or the negative energy) branch at the edge? A natural 
resolution of this puzzle is that the n = Landau level 
is broadened with its bandwidth becoming finite. If so, 
what would be the nature of such bandwidth-broadened 
n = Landau level? 

The quantum mechanical problem of an electron mov- 
ing in lattice under magnetic field is generally known as 
the Azbel-Hofstadter problem named after Azbel |1] , who 
originally proposed the model, and Hofstadter l5], who 
first obtained a numerical solution in the square lattice 
and showed the existence of a self-similar fractal struc- 
ture in energy eigenvalue, dubbed as the Hofstadter but- 
terfly. The actual equations, that need to be solved, 
are known as Harper's equations which are in fact noth- 
ing but the energy eigenvalue equation for the Hamilto- 
nian matrix. By numerically solving Harper's equations, 
the self-similar fractal structure of the Azbel-Hofstadter 
model was found also for various other lattices including 
the triangular and the honeycomb lattice [BHS]. 

In addition to numerical studies solving Harper's equa- 
tion, there have been extensive efforts to obtain analytic 
solutions [TnH26| . The reason for such efforts is multi- 
faceted. For one, many researchers have been curious 
about the very origin of the self-similar fractal structure 
seen in the Hofstadter butterfly and tried to make a con- 
nection to other known systems exhibiting similar fractal 
structures. For another, numerical computations can be 
performed only in the situation where the magnetic flux 
per unit cell, cj), is a rational fraction of the magnetic 
flux quantum, (j}^ = hc/e. Therefore, what happens at 
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irrational fractions can be addressed only by the analytic 
approaches. Perhaps, the most important reason in con- 
nection with experiment is the fact that the numerical 
approach cannot access the weak-field limit where the 
size of the matrix that needs to be diagonalized diverges. 
The weak-field limit is most relevant for experiment since, 
even in the quantum Hall regime, the magnetic flux per 
unit cell is typically much less than 1/100 in units of 
magnetic flux quantum. 

Among various analytic approaches, the Bethe-ansatz 
approach is regarded to be most systematic, where the 
Azbel-Hofstader problem is converted into solving the 
Bethe-ansatz equations whose roots are directly con- 
nected to the energy eigenvalues as well as eigenstates. 
Despite providing such insightful relationship to an in- 
tegrable model, the Bethe-ansatz approach is proven to 
be of little practical use since the Bethe-ansatz equations 
are generally insoluble except for special cases. The use 
of other analytic approaches is also similarly limited. 

In this paper, we develop a method that can be used 
to provide an accurate analytic description of the evolu- 
tion of the n ~ Landau level as a function of magnetic 
field ranging from being arbitrarily weak to moderately 
strong. In this method, it is shown that, for (/)/0o — p/q 
with p and q being coprime positive integers, the central 
band of the Hofstadter butterfly, which is obtained from 
the original 2q x 2q matrix for Harper's equations, is cap- 
tured extremely accurately by diagonalizing the effective 
Hamiltonian matrix with a much reduced size of 2p x 2p 
in the weak-field regime. The central band of the Hofs- 
tadter butterfly is connected with the n ~ Landau level 
in the continuum limit. Actually, this effective Hamilto- 
nian matrix works quite well for (/)/0o as large as 0.3. One 
of the most important discoveries of this work is that, no 
matter how small the magnetic flux per unit cell may 
become, the central Hofstadter band (CHB) always con- 
tains massless Dirac particles whose energy dispersion is 
completely isomorphic to that in the absence of magnetic 
fleld. In fact, by combining the self-similar pattern of 
the central Hofstadter band and some analytic as well as 
numerical results for the zero-energy modes of Harper's 
equations, we conclude that there should be exactly 2q 
Dirac cones in the magnetic Brillouin zone (MBZ) for 
0/00 = p/q with arbitrary p and q. A corollary of this 
result is that there should also be a self-similar occur- 
rence of Landau levels associated with such Dirac cones. 

In order to assess the experimental feasibility of observ- 
ing such massless Dirac particles within the central Hof- 
stadter band, we compute the width of the central Hofs- 
tadter band which, for small 0/0o, is predicted to scale as 
exp(— in units of the energy level spacing between 

the n = and 1 Landau level, y/2hvp/lB- Here, 7 — 
|Cl2(57r/3)|/7r ~ 0.323 and 012(6*) = E^^i sin M)/^^ is 
called the Clausen function. Actually, motivated by an 
intriguing conjecture proposed by Thouless |27j a while 
ago, there has been a long history for addressing how the 
total bandwidth of the Hofstadter butterfly scales as a 
function of magnetic field [3 [T7l[25H5^ . To the best of 



our knowledge, our result is the first report for the scaling 
of the Hofstadter butterfiy bandwidth in the honeycomb 
lattice. Considering difficulties in directly observing the 
Hofstadter butterfiy under magnetic field with typically 
available strength, we believe that a precise measurement 
of the bandwidth itself can be used to infer the existence 
of the Hofstadter butterfly in addition to the Diophantine 
equation for the quantized Hall conductance 33-37]. 

The rest of the paper is organized as follows. In Sec.|llj 
we present the Azbel-Hofstadter model in graphene with 
a particular choice of gauge called the optimal gauge. In 
Sec. |III[ we analyze various properties of the zero-energy 
solutions for Harper's equations, which play a crucial role 
in our effective Hamiltonian method by generating basis 
wave functions for the central Hofstadter band. A pre- 
cise mathematical form of the effective Hamiltonian is 
presented in Sec. |IV[ where it is shown that the resulting 
magnetic band structure provides an excellent agreement 
with that of the central Hofstadter band obtained from 
the original Harper's equations in the weak-fleld regime. 
In Sec. |V] by using such effective Hamiltonian method, 
we carefully analyze the self-similar recursive pattern of 
the central Hofstadter band, which is then combined with 
analytic as well as numerical results for the zero-energy 
modes to show that massless Dirac particles should occur 



under arbitrary magnetic field. We conclude in Sec. VI 



II. AZBEL-HOFSTADTER PROBLEM FOR 
GRAPHENE 

The Azbel-Hofstadter problem is nothing but an en- 
ergy eigenvalue problem of the tight-binding Hamiltonian 
under magnetic field: 
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where tij is the hopping amplitude between the nearest- 
neighboring sites with its phase determined via the 
Peierls substitution, tij = tg — tge'^'^'^^^ , where 0^ = 
f!' A ■ dl and A is the vector potential. Here, to is 
the hopping amplitude in the absence of external mag- 
netic field. For convenience, we now fix the energy scale 
by setting to ~ 1. The physical energy scale can be re- 
stored by re-introducing to, when necessary. While any 
vector potential satisfying the condition that the contour 
integral, ^ A • rfl, around the hexagonal unit cell equals 
the magnetic fiux per unit cell, 0, is legitimate, we take 
a particular choice of the gauge where only one of the 
three 0ij's adjoining the nearest-neighbor carbon pairs is 
set to be non-zero. The situation is illustrated in Fig. [T] 
This gauge is called the optimal gauge since the size of 
the magnetic unit cell (MUC) is optimal with its value 
being qSo for (/)/<t'o = p/Qj where So is the area of a sin- 
gle hexagonal unit cell [HI |37l |38] . Note that the size of 
the magnetic unit cell is doubled in the usual Landau 
gauge [311 inj. 
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FIG. 1: (Color online) Schematic diagram for the gauge used 
in this work. The yellow parallelogram depicts a magnetic 
unit cell (MUC). Magnetic unit cells are denoted by the MUC 
index, a, along the y-direction. Different carbon sites within 
the same magnetic unit cell are distinguished by the dimer in- 
dex, n, and the A/B sublattice index. Note that horizontally- 
connected A and B carbon sites share the same dimer index. 
Red arrows indicate the directions of the paths, along which 
non-zero phases are gained via the Peierls substitution. The 
value of the non-zero Peierls phase is written near each arrow 
while all the other phases are zero. This gauge is called the 
optimal gauge. 



where a denotes the position of a given magnetic unit 
cell along the y-direction and the dimer index, n, indi- 
cates the position of each dimer within the magnetic unit 
cell. The size of the magnetic unit cell is determined by 
the magnetic lattice translation symmetry. For a ratio- 
nal value of the magnetic flux per unit cell in units of 
magnetic flux quantum, 0/0o = p/q^ Harper's equations 
in Eq. ^ and ^ become periodic with respect to the 
diagonal lattice translation operation oi n n + Iq with 
I being an arbitrary integer. Thus, in this situation, the 
magnetic unit cell covers the dimer index ranging from 
uq to no + q — 1 with ng being an arbitrary initial dimer 
index. See Fig. [l] for illustration. 

Harper's equations in Eq. ([2]) and (|3]) can be simpli- 
fied by using the lattice translation symmetry along the 
y-direction. That is to say, the a-dependence can be re- 
moved by defining the crystal momentum, ky, via the 

Bloch theorem, ^pan = V'n^ e*'^«", with ky — ky\/3a. In 
this representation. Harper's equations are given by 

^^V'X = ^«(^.)Vv!-ii. + ^^^.' (4) 

where 

AM = 2e'("'^^''*) cos (nn^^ + 'f^. (6) 



In the optimal guage. Harper's equations can be writ- 
ten as follows: 

Ei^L = + V^L + e'™'^V^?+i,„_i , (2) 

Ei^L = <„+! + V'l + e-^^^("+^)^ V^i,„+i , (3) 



By realizing that the Bloch condition along the diago- 
nal direction, V„fc^ = e''''"'<t>~k,kjn) with <Pk^k^{n) being 
a periodic function of n with period q, is equivalent to 

the boundary condition, "^n+q k ~ ^^'^'^'^'^nk ' ^^'^ "^^'^ 
convert Harper's equations to an eigenvalue problem of 
the following 2q x 2q Hamiltonian matrix: 
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where riQ, the first dimer index for a given magnetic unit not affect the energy eigenvalue. Note that kd is the di- 
cell, can be chosen arbitrarily since the choice of tiq does agonal momentum measured in units of l/^/ia. Figure^ 
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FIG. 2: (Color online) Hofstadter butterfly showing the en- 
ergy eigenvalue, E /to, as a function of the magnetic flux per 
unit cell in units of magnetic flux quantum, 4>/4>o- Here, to 
is the hopping amplitude in the absence of external magnetic 
field. 



shows all energy eigenvalues of the above Hamiltonian 
matrix as a function of the magnetic flux per unit cell, 
in units of magnetic flux quantum, (pQ. This diagram is 
known as the Hofstadter butterfly. Note that our result 
is completely identical to the previous result obtained by 
Rammal using the Landau gauge [7]. 



III. ZERO-ENERGY MODE 

While every energy eigenvalue of the Azbel-Hofstadter 
problem can be in principle found numerically, the size of 
the Hamiltonian matrix, that needs to be diagonalized, 
diverges as 2q x 2q when g — ^ cx) in the weak-field limit of 
(/)/0o = pIi with fixed p. Obviously, a better approach 
is necessary in the weak-field regime. In this work, we 
present an effective Hamiltonian method that can be used 
to provide an accurate analytic description of the central 
band of the Hofstadter butterfly in the weak-field regime. 

Evolving from the n = Landau level, the central Hof- 
stadter band (CHB) is most intriguing since it may hold 
a key to the mysteries of the fractional quantum Hall ef- 
fect (FQHE) in graphene. Note that, while the fractional 
quantum Hall effect has been observed in graphene, its 
detailed properties are not yet fully consistent with cur- 
rent theoretical understanding [HVt54) . For one thing, 
the excitation energy gap, which is the most essential 
physical observable determining the electron transport, 
is orders-of-magnitude smaller than the corresponding 
theoretical predictions. While this discrepancy could 
be explained by various perturbations such as disorder. 
Landau-level mixing, or ripples of the graphene layer, it 
is believed that the conclusive explanation for its true ori- 



gin is still missing. We think that a precise understanding 
of the nature of the central Hofstadter band can serve as 
an important step towards achieving such explanation. 

Our effective Hamiltonian method is based on the ob- 
servation that (i) all energy eigenstates of the central Hof- 
stadter band are well approximated by those of the zero 
energy, which we call the zero-energy modes, and thus (ii) 
a very accurate effective Hamiltonian can be constructed 
by generating basis wave functions from the zero-energy 
modes. In order to facilitate the discussion for how to 
construct the effective Hamiltonian, let us first investi- 
gate various properties of the zero-energy modes in this 
section. Actual construction of the basis wave functions 



is performed in Sec. IV 

For E — Harper's equations in Eq. Q and 
come decoupled between sublattice A and B: 
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(8) 
(9) 



where V'no ^-nd V'no (which are the amplitudes of the 
wave function at n = rip for sublattice A and B, re- 
spectively) can be regarded as simple normalization con- 
stants. Seemingly otherwise, Eq. (|8| and ([9| are not yet 
the solutions for Harper's equations since the momentum 
is not specified. The momentum is fixed by imposing the 
boundary condition, V'n+g = e^'''''^ipn (which is due to the 
Bloch theorem). The situation is a bit unorthodox here 
since the computation is performed in reverse order to the 
conventional scheme where the energy eigenvalue is de- 
termined for a given momenum. In the current scheme, 
we seek for the right momentum corresponding to the 
zero-energy solution. 

To find the right momentum for the zero-energy mode, 
it is convenient to use the following cosine product iden- 
tity: 



n+q 

n 

m— n+l 
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(10) 



where 7pq„ = pn + 1 + {q + l){p — l)/2. The derivation 
of the cosine product identity is given in Appendix 
By using the cosine product identity, one can simplify 



^n+n/^n as foUows: 
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'e^'S" sin ((fey + 7r)|). 
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where Spqn = 2np+ {p+ !)((? + 1) — {q + l)/2. By noting the fohowing: 
that 2np and {p+l){q + l) are always even integers withp 
and q being coprime, one can re- write Eq. ( [TT| ) as follows: 



2e- 



K<?+i)e'^sin((fc, + 7r)|). (12) 



Then, the boundary condition, Vn+q/V'n = e^^'''' , gives 
rise to the following equation for the zero-energy mode 
momentum: 



2 sin I (fc- 



1\ - „»fed9-*^+»f(9+l) 



(13) 



from which ky and can be simultaneously determined. 
First, noting that the magnitude of the left-hand side 
should be unity, one can determine ky by imposing 



3in((fcy + 7r)|) 



(14) 



with j being an integer. The solution of Eq. (14 1, fc*, is 
given by: 



39 
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(15) 
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Then, by inserting Eq. (15) into (13 1, one can determine 
the other momentum for the zero-energy mode, fc^, whose 
value is given as follows: 



kj 



39 
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(16) 
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with / being an integer. Note that, while the preceding 
computation is performed only for sublattice B, it can 
be shown that the zero-energy momentum is exactly the 
same for sublattice A as well. So far, the conclusion is 
that the wave function for the zero-energy mode is pre- 
cisely described by Eq. |8| and ( [9| w ith the appropriate 
momenta given by Eq. (|T5[) and (16 1. 



At this point, it is illuminating to obtain the locations 
of the zero-energy momenta in the magnetic Brillouin 
zone. To this end, let us convert kd in terms of the usual 
cartesian coordinates. Since kd is the diagonal momen- 
tum along which the dimer index, n, increases within a 
given magnetic unit cell index, a, the conversion rule is 
given by 



kd = Viakd = '^ak^ + ^aky, 



(17) 



which, combined with Eq. (15 1 and (16), gives rise to to 



2_n^ 
3q 



(18) 



where I' = 21— j. Figure|3]presents the energy dispersions 
of the central Hofstadter band in the form of contour plot 
for various flux values, where the positions of the zero- 
energy momenta are denoted by little x marks. As one 
can see, the zero-energy momenta occur exactly in the 
same honeycomb pattern as the Dirac points in the ab- 
sence of magnetic field. Actually, it is shown in Sec. IV C 
and |IVD| that, in the weak- field limit, with proper en- 
ergy and momentum re-scaling, the energy dispersion of 
the central Hofstadter band becomes exactly identical to 
that of graphene in the absence of magnetic field, prov- 
ing that the zero-energy modes are, in fact, nothing but 
massless Dirac particles. The energy dispersion remains 
very close to that in the absence of magnetic field even 
when the magnetic flux per unit cell becomes moderately 
large. 

It is interesting to mention that the number of zero- 
energy modes is given by 2q within each magnetic Bril- 
louin zone and this fact is related with the Landau-level 
degeneracy of graphene in the continuum limit. In the 
lattice model, the degeneracy of the Landau level can 
be regarded as the number of different ways of locating 
the wave packet maximum within the magnetic unit cell. 
Since the magnetic unit cell contains 2q carbon atoms, 
the wave packet maximum can have 2q different locations 
and therefore the so-defined Landau-level degeneracy is 
2q, which, in the continuum limit, becomes infinite, or a 
macroscopic number proportional to the system size. 

We now investigate the wave function profile for the 
zero-energy mode. The wave function for the zero-energy 
mode can be computed numerically by using Eq. ([s]) and 
([9]). Figure [4] shows the results for several different flux 
values. One of the most salient features of the exact 
wave function profile is the fact that it is asymmetric 
around its maximum point while, in the continuum limit, 
the zero-energy wave function reduces to the Gaussian 
wave packet (which is the energy eigenstate in the n = 
Landau level) and therefore should be symmetric. As 
one can see from Fig. |4j however, the deviation from the 
Gaussian shape vanishes rather rapidly as the flux per 
unit cell decreases. 

Actually, in the weak- field regime, it is possible to de- 
rive a better analytic approximation for the zero-energy 
wave function than the simple Gaussian. The basic idea 
is, first, to convert the zero-energy wave function repre- 
sented in a product form to a summation form by taking 
the logarithm and, then, to approximate the summation 
with an integral by regarding, a;„ = nTT^/i^o + ky/2, as a 
continuous variable. This procedure is valid when (j)/(f>o is 
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FIG. 3: (Color online) Contour plots for the energy dispersion at various flux values of cj>/cj>o = p/q, with p = 1 and q 
increasing from 1 to 6 between panel (a) through (f). In the figure, the energy dispersions are normalized by their respective 
half bandwidth. The positions of the zero-energy momenta are denoted by little x marks and the magnetic Brillouin zones are 
enclosed by red solid lines. The above energy dispersions are computed by solving either the original Harper's equations or 
the effective Hamiltonian method explained in Sec. |IV| both of which produce essentially the identical results. It is interesting 
to note that the effective Hamiltonian method works well even for p/q — 1 owing to the mirror structure of the Hofstadter 
butterfly, which maps the region near p/q = 1 to the weak-field counterpart. 



small. Relegating the detailed derivation to Appendix [B| 
here, we present the final result: 



\ip„ cx exp 



cx exp 



^, CI2 f 2T:—n + r] 
_2tt(P/4>q \ (j)Q 



(19) 



where rj = ky+7T{(p/(j)o + l) and 012(6*), called the Clausen 
function, is defined such that 012(6*) = J2^=i sin {n9)/in?. 
From now on, let us call the wave function profile given 
by Eq. ([l9]) the Clausen wave packet. As one can see from 
Fig. [4] the Clausen wave packet provides a very accurate 
approximation of the exact results for a wide range of 
flux values. 

To confirm analytically that the Clausen wave packet 
indeed reduces to the Gaussian in the continuum limit, 
it is convenient to use the Landau gauge, in which case 
the Clausen approximation corresponds to the following: 



oc exp 



1 



27r 
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27r — n 

(Po 



(20) 



where k — fc„ 
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TT. Here, we only consider the 



amplitudes in sublattice B since the same process can 
be applied to sublattice A. With the definition of a new 
continuous variable, x = ^(n — 1) (where a is the lattice 



constant), Eq. (20) can be re- written as follows: 



(x) cx exp 



5*0 



CI2 



/2 



-X + K 



(21) 



where 27r0/(/)o = Sq/Iq, with Sq — ^^^a^ being the area 
of the hexagonal unit cell, is used. In the above, k' = 
K + So/ll- 

Now, noting that Eq. (21) is maximized when 



the Clausen function inside the exponential becomes 
minimized, we search for the condition minimizing 
ChiVSax/l^+n'). To this end, it is convenient to use the 
following property of the Clausen function that ChiO) 
has a maximum at 9 = 7r/3 and a minimum at 57r/3 
within a single period between and 2tt. Note that ChiO) 
is a periodic function with 2tt period. Then, one can de- 
termine the maximum position of |?/'^(a;)| as follows: 



1 



/2 



^/3a/2 



fiT: 



(22) 



is the differ- 

3\/3a 

ence between ky and the momentum of one of the two 



where n is an integer and qy = ky 
ence between ky and the momentu 
Dirac points. (Note that, for sublattice A, qy is defined 
as the difference between ky and the momentum of the 
other Dirac point.) Since the Clausen function can be 
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FIG. 4: (Color online) Wave function profiles for the zero- 
energy mode as a function of dimer index, n, at three different 
flux values: (a) 0/0o = p/g = 1/5, (b) 1/20, and (c) 1/100. 
As one can see, at moderate flux values, say, p/q = 1/5 and 
1/20, there is a sizable asymmetry around the maximum po- 
sition. The asymmetry is seen more clearly in contrast to the 
Gaussian wave packet (red dashed lines) which is the exact 
energy eigenstate in the continuum, or weak-field, limit. It is 
important to note that, while the Gaussian wave packet pro- 
vides a poor representation of the exact results (solid lines) at 
moderate flux values, a new analytic expansion method using 
the Clausen function (open circles) works very well for a wide 
range of flux values. 



expanded around its minimum positions as follows: 



Ch{e) - Cl2(57r/3) 



4 ^ 



(23) 



where Cl2(57r/3) = -1.0149 and ^^in - Stt/S + 2^7r 
with /i being an integer, the next step is to expand the 
Clausen wave packet in the vicinity of Xmaxj assuming 
\x — a^maxl/a It is important to note that such 

expansion becomes very accurate when the inverse coef- 
ficient in front of the Clausen function, I'^/Sq, is much 
larger than the deviation of the Clausen function from 
its minimum position. The expansion is given as follows: 




[Ao + X2{x - Xn 

Do 



(24) 



where the linear term vanishes due to the extremum con- 



dition. As shown from the comparison with Eq. (231, the 
zeroth-order coefficient, Aq, is equal to Cl2(57r/3) and the 

second-order coefficient is given by A2 = ^ 



S'o/(2/g), which finally gives rise to the desired result 
that the Clausen wave packet reduces to the usual Gaus- 
sian function of exp {—{x — Xi^^x)^ /2l^). Note that this 
result is exactly the same as the previous result obtained 
by Goerbig and collaborators 40]. 

For later use, it is convenient to compute the max- 
imum as well as the minimum positions of the zero- 
energy wave function for the optimal gauge in the weak- 
field regime. In the case of sublattice B, the maximum 
(minimum) position arises whenever the cosine factor of 

An{ky) in Eq. (6), cos(n7r^ -I- ^) , passes through 1/2 
from above (below) to below (above) as a function of 
dimer index, n. Note that n-T^-^ can be treated roughly 
as a continuous variable so long as 0/(/>o is sufficiently 
small. With the maximum and the minimum position 
denoted as 
follows: 



and nj^ii„, respectively, the result is as 



floor 



= floor 




(25) 



where s is an integer. In the case of sublattice A, it can 
be shown that n^^^ = n^^^ and n^^^^ = n^^^ since the 
cosine factor is multiplied inversely in this case. Finally, 
it is interesting to mention that, in the strong-fleld regime 
where the magnetic flux is in the vicinity of unity, i. e., 
|0/(/)o — 1| <Si 1, the maximum and the minimum-position 
formula is modified as follows: 



. = floor 



nnun,s = floor 



1 




ky 


7r(l - 0/ 


M (3 ^ 


2 


1 


/ 2n 




7r(l - 0/ 


M [ 3 


2 



, (26) 



where s is, again, an integer. 



IV. EFFECTIVE HAMILTONIAN 

In the preceding section, we have carefully investigated 
various aspects of the zero-energy solution for Harper's 
equations. Despite many nice, analytic properties, the 
zero-energy modes alone consist of only a negligible part 
of the entire magnetic Brillouin zone. While all energy 
eigenvalues can be, in principle, computed by solving 
Harper's equations, a brute-force numerical diagonaliza- 
tion is prohibited in the weak-field regime where the size 
of the Hamiltonian matrix quickly diverges. To scan the 
entire Brillouin zone in the weak-field region, it is nec- 
essary to devise a better method. In this section, we 
present such a method using the effective Hamiltonian, 
which provides a very accurate description of the central 
Hofstadter band in the entire Brillouin zone. 
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A. Basis wave functions 

The essence of our effective Hamihonian method hes 
in choosing the right set of basis wave functions most 
relevant to the central Hofstadter band. To do so, it 
is important to note that, for (p/ipo = p/li the central 
Hofstadter band always contains 2p subbands. One way 
of understanding this is, first, to realize that the wave 
function profile of all energy eigenstates comprising the 
central Hofstadter band is more or less identical to that of 
the zero-energy modes in the weak-field limit. Then, from 
Eq. ( [25| , one can see that there should be exactly p local 
maxima for the wave function profile inside the magnetic 
unit cell (while their individual maximum values can be 
different). Now, increases from to 27r 

so that the entire Brillouin zone is covered along the y- 
direction. According to Eq. (251, this process is actually 



identical to decreasing s by unity, which in turn means 
that that the wave function is translated exactly by one 
unit of the distance between the nearest maxima. This 
process covers only 1/p of the whole magnetic unit cell. 
To fill the whole magnetic unit cell, p bands are necessary. 
Since the same is true for both sublattice A and B, there 
should be 2p subbands for the central Hofstadter band. 

Now, we present a scheme for systemically construct- 
ing approximate, but very accurate basis wave functions 
for such 2p subbands. This basis-constructing scheme is 



best explained in the following three steps, (i) First, for a 
given momentum, k = (kx, ky), we compute a trial basis 
wave function by using the zero-energy formula in Eq. Q 
and ([9]). For the time being, let us ignore normalization. 

(ii) We then slice the so-obtained trial wave function into 
equally-spaced p pieces such that each piece contains ex- 
actly one local maximum in the region located between 
two consecutive minima of the trial basis wave function. 
Care must be taken for sublattice A where, according to 
our convention, the boundary of the magnetic unit cell 
sits right on top of one of the wave function maxima and 
thus the piece-wise basis wave function containing such 
maximum is split into two regions separated across the 
magnetic unit cell. In this case, to satisfy the periodic 
boundary condition dictated by the Bloch theorem, we 
multiply an additional phase factor, e*'^'^'^, to the copied 
portion of the wave function amplitude translated from 
the outside to the ending part of the magnetic unit cell. 

(iii) By normalizing the p piece- wise basis wave functions 
separately for each sublattice, we finally obtain 2p basis 
wave functions. Note that the finally obtained basis wave 
functions are orthonormal to each other. See Fig. |5] for 
an illustration of the basis-constructing scheme. 

Explicitly, the basis wave function for sublattice B, 
X^(n) with s ranging from 1 to p, can be written as 
follows: 




for rij^^ _i < n < n^^^^ , 



otherwise. 



(27) 



where cf is the normalization constant. Note that X^(n) 
is the piece-wise basis wave function containing the s-th 
maximum. For sublattice A, the situation is similar ex- 



cept for the special case of s = 1 where the wave function 
maximum is split into two regions across the magnetic 
unit cell: 



J 



„AT-r" 

''I llm=no-|-l 



for no < n < 



for n^ji„ p < n < no + q - 1, 



otherwise, 



(28) 



r 



where ng is the first dimer index in the magnetic unit cell, 
which, according to our convention, is n^^^^^Q + 1. Note 
that the last dimer index is no + <? — 1, which is in turn 



equal to n: 



^^^^^^ ^. In the above, c 
constant. For the other cases wit" 
given similarly to that of sublattice B 



is the normalization 
s =/= 1, the formula is 
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FIG. 5: (Color online) Schematic diagram for the construction of basis wave functions in the case of p/q = 3/152. Basis wave 
functions for the central Hofstadter band can be constructed in three steps: (i) for a given momentum, k = {kx,ky), one 
generates a trial basis wave function according to the zero-energy formula in Eq. ^ and Q, (ii) then, slices the so-obtained 
trial wave function into equally-spaced p pieces such that each piece contains exactly one local maximum, and (iii) finally, 
normalizes the p piece-wise basis wave functions separately for each sublatttice A and B. Refer the text for details. 



1 In 



for 
otherwise, 



(29) 



where, again, is the normalization constant. 

B. Constructing the effective Hamiltonian 

The basic idea behind our effective Hamiltonian 
method is to isolate the Hilbert space near zero energy 
in terms of the basis wave functions constructed in the 
preceding section. With p number of basis wave func- 
tions for each sublattice A and B, say, and with 
fijV — I,-- - ,p, our Hamiltonian can be written as a 
2p X 2p matrix as follows: 

= (jj.Bt ) , (30) 

where H*^ is a p x p matrix whose elements are given 
by 

(H^%. = (x^|H|x^). (31) 

In the above, H is the original Hamiltonian matrix for 
Harper's equations given in Eq. ([7|. Note that all ele- 



ments in the block-diagonal part of H'' are strictly zero 
since H allows only the nearest-neighbor hopping. 



C. Approaching the continuum limit along 

0/00 = 1/g 

The effective Hamiltonian takes the most compact 
form in the case of <^/0o = 1 1 <^■ The reason is that, in this 
case, there is only a single basis wave function for each 
sublattice and thus the size of the effective Hamiltonian 
becomes just 2 x 2 no matter how large q may become. In 
fact, it is important to note that the larger q becomes, the 
more accurate results our effective Hamiltonian method 
provides, as shown later in this section. In addition to 
the mathematical simplicity, the case of 0/0o = 
physically important since taking the large-g limit along 
't>l<i>Q = 1/9 is one of the most natural paths approach- 
ing the continuum limit, via which the central Hofstadter 
band evolves into the n = Landau level. 

With all diagonal elements vanishing (for the reason 
explained in the preceding section), the only non-zero. 
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ofF-diagonal elements of the 2x2 effective Hamiltonian 
are (H^^)ii and its complex conjugate: 

(H^%i = (xf|H|x?) 

— \Xl,no) lAl,rio+^"o6 Xl,no+q-l 
f _ no+q-1 ^ 



} 



(32) 



where C = {Xinn)*XA.n„ ^-^d the cosine product identity 



in Eq. (12 1 is used to obtain the last line. The step con- 
necting between the first and the second line of Eq. ([32]) 
indicates that only a single term from the inner product 
survives. This is due to the fact that all the other terms 
vanish strictly by the very definition of the basis wave 
functions given in Eq. (27), (28), and (29 1, which, in the 
case of 4>l4>o = 1/(7, is simply identical to the zero-energy 
formula in Eq. ([8| and (|9| due to the fact that there is 
only a single maximum in the magnetic unit cell in this 
case. 

Diagonalizing the 2x2 effective Hamiltonian gives rise 
to the following energy eigenvalues. 



£;±^(fc) = ±|c| 



y 1 + 4cos2 (iky ~ fcrf)|) - 4(-l)9 cos ({ky - fed) I) COS ((ky + fed) I) , (33) 
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FIG. 6; (Color online) Overlap integral between the energy 
eigenstates obtained from the effective Hamiltonian method 
and the exact counterparts from the original Harper's equa- 
tions. Stars indicate the averaged value of the overlap integral 
over all crystal momenta within the magnetic Brillouin zone, 
while circles denote individual results for different momenta. 
It is important to note that the overlap integral approaches 
unity very rapidly as 4'l4'o = p/l decreases. 



as well as the corresponding eigenstates. 



1 



(34) 



where 9 is defined such that = 
It is interesting to note that the energy eigenstates are 
always composed of an equal mixture between sublattice 
A and B. 

Figure |6] shows evidence for the validity of the effec- 
tive Hamiltonian method in terms of the overlap inte- 
gral between the eigenstates obtained from the effective 
Hamiltonian and the exact counterparts from the orig- 
inal Harper's equations. As one can see, the overlap is 
very close to unity for all momenta at small flux values 
up to p/q = 0.2. Actually, the overlap is not too bad 
all the way up to p/q = 0.5 when averaged over all crys- 
tal momenta within the magnetic Brillouin zone. Note 
that, for general flux values of p/q, the effective energy 
eigenstates are obtained by solving the 2p x 2p effective 
Hamiltonian. See Sec. lIVDj for details. 

To get m ore physical insight on the energy disper- 
sion in Eq. (33), it is convenient to convert kd in terms 



of the usual cartesian coordinates as done previously in 



Eq. (17). The result is quite illuminating: 



±|CUl -I- 4cos2 (^qaky] -|- 4cos (^qaky) cos {^qak^) 



{q : odd) 



ilCl^l +4cos2 (^^q{aky - ^)) + 4cos (^g(afc, - ^)) cos [fq{ak, - f^)) (g : even) 



, (35) 



which shows that, with proper energy and momentum re- 



scaling, the energy dispersion is, in fact, exactly identical 
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FIG. 7: (Color online) (a) Comparison between the exact 
energy dispersions obtained from the original Harper's equa- 
tions (black dashed lines) and that from the effective Hamil- 
tonian method (red solid line) for various 4i/(j>o = l/?- Note 
that, with proper energy and momentum rescaling, all energy 
dispersions obtained from the effective Hamiltonian at differ- 
ent 4'/4'o ~ 1/9 collapse into a single curve. In the figure, 
the energy dispersions are normalized by their respective half 
bandwidth, W, and the momentum is expressed in units of 
1/qa. The inset shows the path in the magnetic Brillouin 
zone, along which the momentum is scanned. Note that the 
scanning path is chosen such that it passes through the Dirac 
points, (b) Comparison between the exact half bandwidth 
and that obtained from the effective Hamiltonian method as 
a function of <I>/4'Q ~ Vq- 



to that in the absence of magnetic field. Note that, for 
q even, the momentum is shifted by Ak = (1^, ^ ). 

The above energy dispersions were plotted in the form of 
contour graph previously in Fig. [3] for various flux values, 
which shows explicitly that massless Dirac particles exist 
in the central Hofstadter band. 

Figure [7] shows a detailed comparison between the 
exact energy dispersions obtained from the original 
Harper's equations and that from the effective Hamilto- 
nian method for various 0/0o = l/?- It is important to 
note that, with proper energy and momentum re-scaling, 
all energy dispersions obtained from the effective Hamil- 
tonian collapse into a single curve. In the figure, the 
momentum is expressed in units of 1/qa and the energy 



dispersion is normalized by the half bandwidth, W ^ which 
is related with the prefactor, C, via W = 3|C|. As one 
can see from Fig. [T] (a), the agreement between the ex- 
act results for the normalized energy dispersion and that 
from the effective Hamiltonian method is quite good for q 
as small as 3 and becomes perfect quickly as q increases. 
In addition to the re-scaled shape of the energy disper- 
sion, it is shown below that the bandwidth of the energy 
dispersion itself is also captured extremely accurately by 
the effective Hamiltonian method. 

To determine the bandwidth of the energy dispersion, 
it is necessary to compute the prefactor, C, in Eq. (33): 



|C| 



ICaIICe 



(36) 



where the Clausen approximation for the zero-energy 
wave function in Eq. ( 19 1 is used: 



IXl.' 



|CA|e 



CI2 



CI2 27r 



(37) 



Here, Ca and Cb are the normalization constants for sub- 
lattice A and B, respectively. Note that 77 = ky+Tr{(f)/(t)o+ 
1) can be regarded as just a constant for the current pur- 
poses. 

We now need to compute the normalization constants, 
Ca and Cb- First, due to the sublattice symmetry, |Ca| = 
|Cb|, and therefore \C\ — \Cb\^ ■ Mathematically, this is 
a consequence of the property of the Clausen function: 
-Cl2(6') = Cl2(27r - 9). Second, with the substitution 
of = 2TT-$-n + fl, the normalization condition can be 
approximated by the following integral form; 

no+q-l 



|Ce 



ICbP 



n— no 

1 



2tt 



271-0/ 00 Jo 
1 

2^ 



J- 



d0e-^^'^(^'^/^'-A^^ (38) 



where the last line is obtained in the limit of small (f/fpo, 
in which the integrand becomes sharply peaked around 
the minimum position of the Clausen function occurring 
at = 5tt/3 [see Eq. (23)]. In this limit, it is also safe to 
extend the integral range to (— (X),(X)). Following is the 
final result for the half width of the central Hofstadter 
band, W/to: 



W 



1 



-Cl2(5V3) 



(39) 



where we have re-introduced the hopping amplitude, Iq, 
for convenience. Figure [t] (b) shows the comparison be- 
tween the exact half bandwidth and that from the ef- 



fective Hamiltonian method in Eq. ( 39 ) as a function of 
0/00 = l/?! which, as one can see, are in excellent agree- 
ment. It is interesting to note that, in units of the en- 
ergy level spacing between the n = 1 and Landau level. 
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FIG. 8: (Color online) Comparison between the energy dis- 
persions obtained from the original Harper's equations (black 
lines) and those from the effective Hamiltonian method (red 
lines) at various 0/</>o = p/q- For clarity, the energy disper- 
sions from the effective Hamiltonian method are plotted only 
within the window of — 7r/ga < Aky < n/qa. Note that the 
momentum is scanned along the same path as in Fig. [7|(a). 



A = y/2hvp /Ib, the half width of the central Hofstadter 
band becomes simplified as follows: 



W 

'A 



exp 



1 



-Cl2(5^/3) 



(40) 



where Cl2(57r/3) = -1.0149. 



D. General flux 



At general flux, (t>/4>o = p/q, the mathematical expres- 
sion for the energy eigenvalue as well as eigenstate are 
not as simple as those at 0/00 = l/f?, which are given by 



Eq. (33) and (34), respectively, in the preceding section. 



Nevertheless, it is emphasized that, for p/q <C 1, the size 
of the effective Hamiltonian, which is 2p x 2p, is much re- 
duced from that of the original Harper's equation, which 
is 2q X 2q. This means that the fine self-similar structures 
of the central Hofstadter band in the weak-field regime 
can be computed in a much efficient manner. As shown 
in the following section, this, combined with some ana- 
lytic results obtained at 0/0o — l/q, in turn enables us 
to make a prediction that massless Dirac particles should 
occur under arbitrary magnetic field. 

Postponing the detailed discussion to the following sec- 
tion, here, we present the comparison between the re- 
sults obtained from the effective Hamiltonian method 
and those from the original Harper's equations for gen- 
eral 0/</>o =p/q- Figure [s] provides numerical results for 
the energy dispersion at various flux values in compar- 
ison with those from the effective Hamiltonian method. 
As one can see, the agreement is excellent not only for 




FIG. 9: (Color online) A sequence of zoomed views for the 
Hofstadter butterfly in graphene showing various self-similar 
recursive patterns. Note that a fan of narrow energy bands 
are emanated from each single-band boundary flux (SBF), 
0SBF, which, as indicated by blue guiding curves, scale as 
sgn(n)^y\n{(f) — 0sbf)| with n being an integer. This scaling 
behavior is a signature of the formation of recursive Landau 
levels associated with self-similarly occurring massless Dirac 
particles. 



the bands near zero energy, but also for the entire 2p 
bands within the central Hofstadter band. 



V. SELF-SIMILAR OCCURRENCE OF 
MASSLESS DIRAC PARTICLES 

It is mentioned in the preceding section that the ef- 
fective Hamiltonian method can help reveal the fine self- 
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similar structures of the central Hofstadter band in the 
weak-field regime much efficiently. The results obtained 
from the effective Hamiltonian method are shown in 
Fig. [9j which provides a sequence of zoomed views un- 
veiling the self-similar recursive patterns of the central 
Hofstadter band. 

One of the most salient features of the Hofstadter but- 
terfly seen in Fig. |9] is that the central Hofstadter band 
is partitioned by a series of special flux values, at which 
the central Hofstadter band is entirely composed of a sin- 
gle band appearing as a vertical line in the figure. Note 
that, for example, in the top-most panel of Fig. |9] the 
central Hofstadter band is partitioned by a series of ver- 
tical lines occurring at 4>/4'a = 1/9 ^-nd 1 — \/q with q 
being a positive integer. Similar patterns are observed 
in subsequently zoomed views. For later convenience, we 
name the flux belonging to this series of special flux val- 
ues as the single-band boundary flux (SBF). 

In fact, owing to the self-similar recursive pattern of 
the central Hofstadter band, it is convenient to coin the 
name the "n-th level" central Hofstadter band (CHB) 
and the corresponding "n-th level" single-band bound- 
ary flux (SBF). In this nomenclature, the first level SBF 
values are given by 1 /q and \ — l/q with q being a posi- 
tive integer. Meanwhile, the second panel of Fig. |9] shows 
that the second level SBF values are given by 5/49, 4/39, 
3/29, 2/19, 3/28, 4/37, 5/46, and so on. The third level 
SBF values can be determined similarly from the third 
panel. At this point, it is natural to ask the question if 
there is a rule for the SBF values and, if so, what math- 
ematical form it takes. 

The answer is that, indeed, there is a single rule for all 
SBF values, the mathematical form of which is given as 
follows: 



7^l -I- S2 



{-If 



(41) 



n2 + s-i + 



7^3 H 



where n^, a positive integer larger than 2, and Si, either 
or 1, are determined by the following recursion rule. 
Suppose that / is one of the SBF values. Then, we first 
define /o = /. If floor(l//o) > 2, we set ni = floor(l//o) 
and Si = 0. Otherwise, i. e., if floor(l//o) = 1, we set 
ni = floor[l/(l — /o)] and si = 1. As the next recursion 
step, we then define /i = l//o — «i for the former and 
1/(1 — /o) — ni for the latter case. We now repeat the 
same procedure to determine n2 and S2 from /i. This 
procedure can be continued until we get fn = with n 
indicating that / is the n-th level SBF. 

It is instructive to explain the above rule by using an 
example. As an example, let us take '/>/0o — 4/39, which 
is one of the second-level SBF values. According to the 
above-mentioned rule, we first define /o = 4/39. Since 
l//o = 39/4 = 9-1-3/4, Til = 9, si = 0, and subsequently 
/i = 3/4. Now that l//i =4/3 = 1 + 1/3, we have to 
set 712 = floor[l/(l — /i)] ~ 4, in which case S2 = 1. The 



recursion steps terminate at the second level since /2 = 0. 
In conclusion, / = 4/39 can be expressed as follows: 



/ = 4/39 



1 



(-1) 
4 



(42) 



It is now convenient to devise a simplified notation 
scheme where the SBF is represented by a sequence of 
Hi along with whether is or 1. One way of denoting 
the fact that = 1 is to put a bar on top of the corre- 
sponding rii. In this notation, / = 4/39 = (9,4). Similar 
computations can be performed to show that / = 19/186 
and 17/166, which are among the third-level SBF values 
shown in the third panel in Fig. [9j are represented by 
(9,4,4) and (9,4,4), respectively. On the other hand, 
/ = 91/891, which is one of the fourth-level SBF values 
shown in the fourth panel in Fig.|9j is given by (9, 4, 4, 4). 

By knowing the continued-fraction representation of a 
given SBF value, /, one can extract two important pieces 
of information. First, how many n.^'s exist indicates the 
level of / as a SBF value. Second, more importantly, 
provided that / is the m-th level SBF, / is related to 
the first- level SBF occurring at l/rim (or 1 — l/rim via 
the reflection symmetry). For example, / = 19/186 = 
(9, 4, 4) has four n^'s and the last integer is 4, which tells 
us that / = 19/186 is the fourth-level SBF related to the 
first-level SBF occurring 1/4. 

Once the relationship between a given SBF and its 
first-level counterpart is established, there is a far- 
reaching consequence. To understand this, it is impor- 
tant to note that (i) the first-level SBF values are al- 
ways either l/q or I — l/q with q being a positive in- 
teger and (ii) for 4'/(j)a — l/q and 1 — l/q, the energy 
dispersion is isomorphic to that in the absence of mag- 
netic field, as proven in Sec. IV C Therefore, if all SBF 



values are related to their respective first-level counter- 
parts, the energy dispersion at all SBF values should also 
be isomorphic to that in the absence of magnetic field. 
In other words, massless Dirac particles should exist at 
all SBF values. In fact, since all rational fractions can be 
represented by a continued fraction via Eq. 



(411 



less Dirac particles should exist at all rational flux val- 
ues. This conclusion is supported by explicit numerical 
results obtained from both the original Harper's equa- 
tions and the effective Hamiltonian method, which show 
that the energy dispersion is indeed isomorphic to that 
of graphene in the absence of magnetic field. This is, 
also, fully consistent with an analytic result that zero- 
energy modes always exist for general (t>/4>o = p/q as 
shown in Sec. 



HI 



Moreover, since any irrational number 
can be represented as a continued fraction with an infi- 
nite number of levels, the energy dispersion at irrational 
flux values can be regarded as that of massless Dirac par- 
ticles in the limit where the energy scale goes to zero. In 
this sense, we arrive at the final conclusion that, however 
small their energy scale may be, massless Dirac particles 
should exist at all flux values, rational or irrational. 
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A corollary of the above conclusion is that the cen- 
tral Hofstadter band should also contain a self-similar 
structure of recursive Landau levels associated with those 
self-similarly occurring massless Dirac particles. Figure[9] 
shows that each single-band boundary flux (SBF), (/)sbf, 
indeed emanates a fan of narrow energy bands which, 
as indicated by blue guiding curves in the figure, scale 
as sgn{n) \n{(j) — (psBp)] with n being an integer. This 
scaling behavior is a signature of the formation of recur- 
sive Landau levels. 



VI. CONCLUSION 

In this paper, we develop an effective Hamiltonian 
method that can be used to provide an accurate analytic 
description of the central Hofstadter band in graphene 
much more efficiently than directly solving the original 
Harper's equations in the weak-field regime. The source 
of the efliciency is due to the fact that, in the weak-field 
regime where the magnetic flux per unit cell in units of 
magnetic flux quantum, 0/00 = p/q <C 1, the size of the 
effective Hamiltonian is given by 2p x 2p, which is greatly 
reduced from that of the original Hamiltonian, 2q x 2q. 
The benefit of using the effective Hamiltonian method is 
maximized at 0/00 = l/?: where the size of the effective 
Hamiltonian remains to be 2 x 2 no matter how large q 
may become. Actually, the advantage of using the effec- 
tive Hamiltonian is not simply due to the reduction of the 
matrix size, but rather the separation of the low-energy 
sector. It is important to note that solving the original 
Harper's equations generates unreliable, noisy data be- 
low certain small flux values where the low-energy sector 
becomes so narrow that the energy resolution falls below 
numerical accuracy. 

By using such effective Hamiltonian method, we show 
explicitly that the energy dispersion is isomorphic to that 
in the absence of magnetic field for all flux values satis- 
fying 0/00 = l/q, which in turn indicates that massless 
Dirac particles should exist no matter how small the mag- 
netic flux may become. In fact, by combing numerical 
results showing the self-similar recursive structure of the 
central Hofstadter band, we conclude that massless Dirac 
particles should occur under arbitrary magnetic flux. If 
so, as a corollary, the central Hofstadter band should also 
contain a self-similar structure of recursive Landau levels. 

As a useful by-product of the effective Hamiltonian 
method, we are also able to compute the width of the 
central Hofstadter band as a function of magnetic field, 
which can be used to assess the experimental feasibility of 
actually observing massless Dirac particles inside the cen- 
tral Hofstadter band. In units of the energy level spacing 
between the n = 1 and Landau level, A = ^/2Hvf/Ib, 
where vp the Fermi velocity at Dirac point and is 
the magnetic length, we show that the width of the cen- 
tral Hofstadter band is given by W/A — exp (—7^) 
with 7 = |Cl2(57r/3)|/7r ~ 0.323. 

Finally, we mention that the above effective Hamilto- 



nian method is not applicable in the square lattice. The 
reason is as follows. The validity of the effective Hamil- 
tonian method depends crucially on the fact that the 
zero-energy wave function has a well localized shape with 
exponentially negligible tails so that it can be safely split 
into linearly independent pieces with each forming the 
basis wave functions for the effective Hamiltonian. No 
such simplification is possible in the square lattice where 
the zero-energy wave functions are extended all over the 
magnetic unit cell. The situation is not improved in the 
case of non-zero energy states, whose wave function forms 
are no longer given by a simple product form and thus 
prohibit a systematic construction of the analytic basis 
wave functions from the outset. 
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Appendix A: Derivation of the cosine product 
identity 

In this section of Appendix, we prove the following 
cosine product identity: 



n 

rn—n+l 



cos I TOTT^ + aj = sm \^\^a + - j q 



(Al) 



where 7pq„ ^ pn + 1 + (q + l){p — l)/2. Here, p and q 
are coprime natural numbers. 



We begin by multiplying the both sides of Eq. (Al) 
with 2*^, in which case the left-hand side becomes 



n+q 

n 

m—n+l 
n+q 



2 COS I rrnr — ha 

q 



m—n+l 

n+g n-\-q 



m—n+l 



_ 'if p(2n+iy+l)— zag 



m— Ti+1 



n+Q 
m—n+l 



(A2) 



N ow, let us consider the product in the last line of 
Eq. Jli, 5 = n::i'„+i [l + e<2™-t+2")l, whose loga- 
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rithm is written as follows: 

n+q 



m— n+l 

n+q oo / 
_ I -'-J ^i(2m7r£+2a)s 



(A3) 



m— n+l s— 1 



where the Taylor expansion of the logarithm, In (1 + a;) = 

"^"^^i — used. Note that the above Taylor 
expansion of the logarithm is valid for |a::| < 1 with 
exception of x = —1. This condition is satisfied for 
X = e<2™^f +2a) ^j^^ggg g»(2m^E+2a) ^ Fortunately, 

in the case when there is such m that e'^^™'^?"''^") = 



1, Eq. (Al) is automatically satisfied with the both 



sides becoming simultaneously zero. The reason is that 

gi(2m,rE+2a) ^ _^ ^niTT^ + 2a = {21 - 1)tT with I 

being an integer, which is in turn equivalent to 



(a+ |] f? = {ql-pm)n. 



(A4) 



Since ql—pm is an integer, the right-hand side of Eq. ( Al ) 



becomes zero. It is shown in Eq. (A2) that the left-hand 
side also vanishes when g»(2™7r^+2a) _ _ Therefore, it 
can be concluded that the Taylor expansion in the above 
can be safely used. 



Then, Eq. (A3) can be simplified as follows: 

n+q 



-e 2^ e 

ni—n+l 

\.s-l 



s = l 



s^O (mod q) 

i: ^ 

s—0 (mod q) 



z2m7r - s 



{-ly 



i2cts 



1 



1 



e i27r^s(n+l) 
227r — s 



E 
E 



-1)"'-! 

~l 



^i2aql 



E 



(-1) 



(-1 



In 



(A5) 



where the last step in the above equation is obtained 
when e*^["«+7(«+^)] ^ -1. Fortunately, this condition is 
identical to the previous one that there is no s uch m sat- 
isfying e'(^"''t+^") = -1 as described in Eq. (A4). The 
reason is as follows. First, e*^["''+5 •^'+-^^1 ^ -1 indicates 



that (a -|- Tr/2)q = kn with k being an integer. Now that 
p and q are coprime, there should exist integers, ni and 
ri2, such that nip + 7i2<z — 1 according to Bezout's iden- 
tity, which means, in turn, that any integer, say fc, can be 
re-written as {kni)p + [kn2)q- The comparison between 
this condition and that in Eq. (A4) shows that they are 



in fact identical since one can always choose / = kn2 and 
m = —kni. 



Exponentiating the both sides of Eq. ( A5 ) gives rise to 
the following result: 

5^1+ gi2[Qg+f (g+l)] 



cos i^aq + -(q + l) 



(A6) 



By using this result, one can then show that Eq. (A2) 
becomes as follows: 



n+q 

n 

m— n+l 



2 cos I mTT- 

q 



2g-™[pn+l + (q+l)(p-l)/2] gjj. 



::^2e"^'"'" sin(^(^c 



(A7) 



where 7pg„ = pn + l + {q + l){p—l)/2. Dividi ng th e both 
sides of Eq. (IX7| by 2"? finally resuhs in Eq. (IaT). 



Appendix B: Clausen approximation for the 
zero-energy mode 

In this section of Appendix, we derive the analytic ex- 
pression for the wave function profile of the zero-energy 
mode, which becomes exact in the weak-field limit, and 
provides a very approximation to the exact solution at 
moderately small flux values. For completeness, here, we 
consider both the optimal and the Landau gauge. 

In the case of the optimal gauge, let us begin with the 
following Harper's equation for the zero-energy mode in 
sublattice B: 



n 

n [-^M 



(Bl) 



coslTOTT^-l-^). Tak- 



where Am{ky) = 2e V ''J cos 

ing the absolute value and the logarithm of the both sides 
of Eq. (Bl ) gives rise to the following: 



In 



El- 



cos TOTT- 



(f) ky 



(B2) 



In the weak-field limit when (j)/(j)o <C 1, one can approx- 



imate the summation in the right-hand side of Eq. ( B2 ) 
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with an integral via the substitution oi x — rmr4^ + 

and dx = t^-^- That is to say, by using the midpoint 
rectangle method, one can approximate the above sum- 
mation as follows: 



In 



1 


/"■ 


+ 


2 00 


7r0 




■2 


A. 

00 


1 


r 


+ 


IT ^ 

2 00 






2" 


00 




r 


+ 


2 00 


7r0 




2" 


00 



In (2 1 cos a; I) 



rfa; In 1 + e 



oo 



1 sin [s{2x + tt)] 

^ s=l 



_Y cos(2sa;) 
s 

2 00 



2s2 



1 

"27rfi 



CI2 ( 22;„ + TT-^ + TT 



CI2 2X1 - TT^- +7r 



(B3) 



Now, let us switch gears to the Landau gauge, A = 
(0,Bx). In the Landau gauge, the hopping amplitude 
gains the following phase whose value is determined by 
the line integral between the nearest neighboring sites. 



[2(a'-a)-(-in(5„„.0„ 
[2(a'-a) + (-l)"]5„„,0„ (B6) 



where 0„ = ^ (ri/2 — 5/12). As before, n is the dimer 
index and a labels a unit cell along the y-direction (See 
Fig.[l]). 

At this point, it is convenient to consider a semi-infinite 
configuration of graphene with a zigzag edge, in which 
case the wave function amplitude on one of the sublat- 
tices can be chosen to be identically zero. Defining sub- 
lattice B as the one with non-zero wave function ampli- 
tudes, one can show that the wave function amplitude in 
sublattice B is given as follows: 



where x„ = n-^^'^ + ^ and xi = Xn=i- Note that 
012(0) — X^i^Li sin (n6')/n^ is called the Clausen func- 
tion. Neglecting the proportionality constant which is 
independent of n, we arrive at the final result: 



oc exp 



(B4) 



where rj — ky + 7r((/</(/)o + 1)- By noting that Harper's 
equation for sublattice A is simply the inverse of that for 
sublattice B, one can obtain the following expression for 
the wave function profile in sublattice A: 



V'„ oc exp 



1 



27r 



-CI9 



2TT — n 

90 



(B5) 



n 



-2 cos mn 



(j) ky 



(B7) 



Since the above formula is basically identical to that of 
the optimal gauge in Eq. (Bl|, the same computation 



procedure previously applied in the optimal gauge can 
be performed to show that, in the weak-field limit. 



\^n+i \ « exp 



where k = ky — ^ 



, — CI2 2n—n + K 



(B8) 
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